/*

  Simple cross-product, outputs in correct triangular form.

  - Memory reads are float4s, half cooalescing
  - No memory coalescing in writes
  - All threads read in all data required
  - Each thread works on a 2x2 square of data

 */

CUBE_KERNEL(simple2x1, cuFloatComplex *array_d, cuFloatComplex *product_d, 
	    const int write)
{
  CUBE_START;

  //set frequency number from blockIdx.y
  int f = blockIdx.y;

  int Row, Col, blockX, blockY;
  CUBE_DEVICE_CALL(findPosition, Col, Row, blockX, blockY);
  Row*=2;

  //instantiate sum variables
  float sum11XXreal = 0.0, sum11XXimag = 0.0;
  float sum11XYreal = 0.0, sum11XYimag = 0.0;
  float sum11YXreal = 0.0, sum11YXimag = 0.0;
  float sum11YYreal = 0.0, sum11YYimag = 0.0;
  float sum21XXreal = 0.0, sum21XXimag = 0.0;
  float sum21XYreal = 0.0, sum21XYimag = 0.0;
  float sum21YXreal = 0.0, sum21YXimag = 0.0;
  float sum21YYreal = 0.0, sum21YYimag = 0.0;

  //  #pragma unroll 50
  for(int t=0; t<NTIME; t++){
    float4 col1 = ((float4*)array_d)[t*NFREQUENCY*NSTATION + f*NSTATION + Col];
    float4 row1 = ((float4*)array_d)[t*NFREQUENCY*NSTATION + f*NSTATION + Row];
    float4 row2 = ((float4*)array_d)[t*NFREQUENCY*NSTATION + f*NSTATION + Row + 1];
    CUBE_ADD_BYTES(48);

    sum11XXreal += row1.x * col1.x;
    sum11XXreal += row1.y * col1.y;
    sum11XXimag += row1.y * col1.x;
    sum11XXimag -= row1.x * col1.y;
    sum11XYreal += row1.x * col1.z;
    sum11XYreal += row1.y * col1.w;
    sum11XYimag += row1.y * col1.z;
    sum11XYimag -= row1.x * col1.w;
    sum11YXreal += row1.z * col1.x;
    sum11YXreal += row1.w * col1.y;
    sum11YXimag += row1.w * col1.x;
    sum11YXimag -= row1.z * col1.y;
    sum11YYreal += row1.z * col1.z;
    sum11YYreal += row1.w * col1.w;
    sum11YYimag += row1.w * col1.z;
    sum11YYimag -= row1.z * col1.w;
    sum21XXreal += row2.x * col1.x;
    sum21XXreal += row2.y * col1.y;
    sum21XXimag += row2.y * col1.x;
    sum21XXimag -= row2.x * col1.y;
    sum21XYreal += row2.x * col1.z;
    sum21XYreal += row2.y * col1.w;
    sum21XYimag += row2.y * col1.z;
    sum21XYimag -= row2.x * col1.w;
    sum21YXreal += row2.z * col1.x;
    sum21YXreal += row2.w * col1.y;
    sum21YXimag += row2.w * col1.x;
    sum21YXimag -= row2.z * col1.y;
    sum21YYreal += row2.z * col1.z;
    sum21YYreal += row2.w * col1.w;
    sum21YYimag += row2.w * col1.z;
    sum21YYimag -= row2.z * col1.w;

  /*
    sum11XXreal += row1.x * col1.x + row1.y * col1.y;
    sum11XXimag += row1.y * col1.x - row1.x * col1.y;
    sum11XYreal += row1.x * col1.z + row1.y * col1.w;
    sum11XYimag += row1.y * col1.z - row1.x * col1.w;
    sum11YXreal += row1.z * col1.x + row1.w * col1.y;
    sum11YXimag += row1.w * col1.x - row1.z * col1.y;
    sum11YYreal += row1.z * col1.z + row1.w * col1.w;
    sum11YYimag += row1.w * col1.z - row1.z * col1.w;

    sum21XXreal += row2.x * col1.x + row2.y * col1.y;
    sum21XXimag += row2.y * col1.x - row2.x * col1.y;
    sum21XYreal += row2.x * col1.z + row2.y * col1.w;
    sum21XYimag += row2.y * col1.z - row2.x * col1.w;
    sum21YXreal += row2.z * col1.x + row2.w * col1.y;
    sum21YXimag += row2.w * col1.x - row2.z * col1.y;
    sum21YYreal += row2.z * col1.z + row2.w * col1.w;
    sum21YYimag += row2.w * col1.z - row2.z * col1.w; */
  } 

  if (write)
    CUBE_DEVICE_CALL(write2x1, Col, Row, product_d, sum11XXreal, sum11XXimag, sum11XYreal, 
		     sum11XYimag, sum11YXreal, sum11YXimag, sum11YYreal, sum11YYimag, 
		     sum21XXreal, sum21XXimag, sum21XYreal, sum21XYimag, sum21YXreal, 
		     sum21YXimag, sum21YYreal, sum21YYimag);

  if (Col < Row+1) {
    CUBE_ADD_BYTES(64);
    CUBE_ADD_FLOPS(NTIME*64);
  } else {
    CUBE_ADD_BYTES(32);
    CUBE_ADD_FLOPS(NTIME*32);
  }

  CUBE_END;
}
